C.................... RVN2PB.FOR;1 ..........22-JUL-1991 09:42:48.04 
      SUBROUTINE VIBN2P(JP,ID,ISW,N2,O,NE,OP2P,OP2D,RTS
     >  ,QX2S,QA2P,QB2S,O2,N2TOT,ZR)
C... THIS PROGRAM WAS ORIGINALLY WRITTEN BY N ORSINI.  FOR DETAILS SEE HIS
C... PH.D. THESIS 1977(UNIVERSITY OF MICHIGAN). PAGE NUMBERS IN THIS PROGRAM
C... REFER TO THESIS PAGE NUMBERS.
C... The switch ISW is for printing the vibrational populations off N2+
C... and switching the destination of O+(2D) from the X to the A state.
C... ISW=+1 : O+(2D) -> X ; print populations
C... ISW=-1 : O+(2D) -> A ; "        "
C... ISW >+1  O+(2D) -> X ; ISW < -1 (or 0) O+(2D) -> A but no print
C******** In this file, O+(2D) -> A now
C
      IMPLICIT REAL (A-H,N-Z)
      DOUBLE PRECISION N2,O,NE,OP2P,OP2D,RTS(199),O2,N2TOT
     > ,ZR,QX2S,QA2P,QB2S
      DOUBLE PRECISION AWORK(51,51),BWORK(51),DC(51),RK(9),CC(15)
      DOUBLE PRECISION RL(9),PNOMET,ASTOR(15,15),R42,RAV,PROD,
     >  ALOSS,CHN2P,ALR(22),ASOL(15),EMRAT(15),TVBN2P,TVSAV,CEX(9),
     >  TV,PR(22)
      REAL KQ
      DIMENSION A(5,5,4),GG(5,5,4),Q(5,5,4),GH(4),VFAC(5),GTV(5)

      DATA GTV,VFAC,GH/3373.,3373.,6664.,9934.,13162., 5*50.0,  4*1.0/
      DATA    KQ  , TV  ,IMAX,JTI,PFAC,TVBN2P
     > /  5.0E-10 ,1000., 14 , 0 , 1.0,3000.0 /
      DATA BWORK/51*0.0D0/,AWORK/2601*0.0D0/
C
C.........  THE A'S,GG'S AND Q'S COME FROM C.1-C.4 PAGES 143-148
       DATA A/50*0,9.53E4,1.011E5,6.05E4,2.75E4,1.08E4,3.34E4,
     > 7.82E3,6.93E4,8.61E4,6.07E4,3.5E3,2.82E4,2.09E3,2.53E4,
     > 7.31E4,1.0E2,6.1E3,1.5E4,1.34E4,3.03E3,0.0,2.5E2,6.78E3,5.14E3,
     > 1.99E4,1.1E7,6.7E6,1.27E6,0.0,0.0,3.35E6,3.87E6,9.12E6,2.98E6,
     > 0.0,6.75E5,3.81E6,8.99E5,9.46E6,0.0,1.11E5,1.32E6,3.13E6,3.8E4,
     > 0.0,1.61E4,3.12E6,1.71E6,2.18E6,0.0/
C
       DATA GG/50*0,6.42E-2,5.2E-2,1.18E-2,5.13E-4,0.0E-1,4.02E-2,
     >  5.7E-3,8.4E-2,2.5E-2,1.7E-3,1.5E-2,3.0E-2,1.6E-3,2.6E-2,
     >  3.0E-2,8.31E-3,2.2E-2,1.1E-2,1.1E-2,1.1E-2,1.1E-3,2.0E-2,
     >  2.1E-2,1.6E-3,1.9E-2,7.26E-2,6.25E-2,2.66E-2,8.05E-3,0.0E-1,
     >  2.48E-2,2.81E-2,7.34E-2,9.56E-2,1.98E-2,4.1E-3,3.38E-2,8.1E-3,
     >  6.2E-2,6.1E-2,0.0E-1,9.1E-3,3.6E-2,7*0.0/
C
       DATA Q/25*0,9.0236E-1,9.4136E-2,3.4381E-3,6.8513E-5,8.2420E-7
     > ,9.0636E-2,7.2014E-1,1.7826E-1,1.0666E-2,3.0295E-4
     > ,6.5080E-3,1.6493E-1,5.5520E-1,2.5061E-1,2.1923E-2
     > ,4.5437E-4,1.8798E-2,2.2214E-1,4.0980E-1,3.0965E-1
     > ,3.5005E-5,1.8024E-3,3.5913E-2,2.6203E-1,2.8570E-1
     > ,2.4468E-1,3.7227E-1,2.5160E-1,9.9998E-2,2.6080E-2
     > ,3.1069E-1,4.1375E-2,8.1429E-2,2.5650E-1,2.0335E-1
     > ,2.2530E-1,3.5985E-2,1.6578E-1,2.6916E-3,1.2426E-1
     > ,1.2358E-1,1.4331E-1,1.8940E-2,1.1328E-1,7.4351E-2
     > ,5.7340E-2,1.6185E-1,2.2750E-2,9.5849E-2,2.2037E-2
     > ,8.9119E-1,9.7304E-2,1.0474E-2,9.4914E-4,7.2849E-5
     > ,1.0703E-1,7.1633E-1,1.4924E-1,2.4218E-2,2.8875E-3
     > ,1.7514E-3,1.8270E-1,5.9899E-1,1.7279E-1,3.7606E-2
     > ,2.6855E-5,3.5203E-3,2.3645E-1,5.2279E-1,1.7860E-1,5*0/
C
C....... If production is small return
      IF(QX2S.LE.1.0E-10) THEN
         DO 58 J=1,IMAX
 58      BWORK(J)=0.0
         TVBN2P=0.0
         IF(IABS(ISW).EQ.1) WRITE(ID,299) NINT(ZR),(BWORK(J),J=1,IMAX)
     >      ,TVBN2P
         RETURN
      ENDIF
C........ Switch O+(2D) from X to A state. ******* always on now **
C.....      IF(ISW.LE.0) THEN
      PFAC=0.0
      DO 5 I=1,5
 5    VFAC(I)=1.0
C.....      ENDIF
 
      IMAXP1=IMAX+1
      JTI=JTI+1
C     JKL=JP+2
C     JW=((JP+2)/4)*4
C
       IF(JTI.NE.1) GO TO 50
C.......... NEXT LINE PREVENTS OUTPUT, DELETE FOR OUTPUT.........
      IF(ISW.NE.-99) GO TO 50
       DO 30 J=1,5
 30    WRITE(ID,90)(A(I,J,3),I=1,5),(A(I,J,4),I=1,5)
       DO 40 J=1,5
 40    WRITE(ID,90)(GG(J,I,3),I=1,5),(GG(J,I,4),I=1,5)
       DO 35 K=2,4
       DO 35 J=1,5
 35    WRITE(ID,90) (Q(I,J,K),I=1,5)
 50    CONTINUE
C
C........ TVBN2P= vibrational temperature of N2+
      TVBN2P=300.
      ITV=0
C........... loop on TV
C 23  CONTINUE
      TVSAV=TVBN2P
C
      DO 900 I=1,IMAXP1
      DC(I)=0.0D0
      CC(I)=0.0D0
      ASOL(I)=0.0D0
      DO 900 J=1,IMAXP1
  900 AWORK(I,J)=0.0D0
C
C........ RK=RTS(11) = N2+ + e recombination rate (Abdou). RL = N2+ + O
C........ rate. TVBN2P= vibrational temperature of N2+
C........ Rates may be adjusted for different vibrational levels
      RK(1)=RTS(11)
      RL(1)=RTS(42)*RTS(10)
      DO 92 I=1,4
      RL(I+1)=GH(I)*RL(1)*VFAC(I)
 92   RK(I+1)=RK(1) * 1.0
C
C... DIAGONAL COEFFICIENT OF THE MATRIX ON PAGE 149.GG'S AND
C...  A'S ARE EXCITATION AND DEEXCITATION LOSSES
      DO 170 J=1,5
      DO 170 I=1,5
      DC(J)=DC(J)+GG(J,I,3)+GG(J,I,4)
      DC(J+5)=DC(J+5)+A(J,I,3)
 170  DC(J+10)=DC(J+10)+A(J,I,4)
C
C... COEFFICIENTS OF MATRIX ON PAGE 149.  CHEMICAL LOSSES ARE ADDED TO N2*DC'S
C...  N2+ + O -->NO+ + N [RTS(10)]:N2+ + O --> N2 + O+(4S) [RL(I),RTS(42)].
      DO 180 I=1,5
      DC(I)=DC(I)+O*(RTS(10)+RL(I))+RK(I)*NE+RTS(17)*O2
C... DEACTIVATION OF N2+(X2S)V=1,4 BY COLLISION WITH N2
      AWORK(1,I)=KQ*N2
      IF(I.GT.1) DC(I)=DC(I)+KQ*N2
      DO 180 J=1,5
      AWORK(J,I+5)=A(I,J,3)
      AWORK(J,I+10)=A(I,J,4)
      AWORK(J+5,I)=GG(I,J,3)
 180  AWORK(J+10,I)=GG(I,J,4)
C
C... FILL IN THE DIAGONAL ELEMENTS OF MATRIX PAGE 149
      DO 190 I=1,IMAX
 190  AWORK(I,I)=-DC(I)
C
C... CALCULATION OF C'S ON PAGE 141
      CEX(1)=(1-EXP(-GTV(1)/TV))
      DO 195 I=2,5
 195  CEX(I)=CEX(1)*(EXP(-GTV(I)/TV))
C
C... CALCULATION OF Q X C ON PAGE 142
      DO 220 I=1,5
      DO 220 J=1,5
      CC(I)=CC(I)+Q(J,I,2)*CEX(J)
      CC(I+5)=CC(I+5)+Q(J,I,3)*CEX(J)
 220  CC(I+10)=CC(I+10)+Q(J,I,4)*CEX(J)
C
C........  PNOMET= PROD OF N2+ FROM META O.....PFAC=1 PUTS PNOMET INTO
C........  N2+(X2S),PFAC=0 PUTS IT INTO N2+(A2P)
C
      PNOMET=(RTS(19)*OP2D+RTS(20)*OP2P)*N2
C... RIGHT HAND SIDE OF MATRIX ON PAGE 149
      DO 230 J=1,5
      AWORK(J,15)=-(QX2S+PFAC*PNOMET)*CC(J)
      AWORK(J+5,15)=-QA2P*CC(J+5)
 230  AWORK(J+10,15)=-QB2S*CC(J+10)
C.....
      AWORK(6,15)=AWORK(6,15)+0.5*(PFAC-1)*PNOMET
      AWORK(7,15)=AWORK(7,15)+0.5*(PFAC-1)*PNOMET
C
C...  IF(ISW.NE.0.AND.JW.EQ.JKL) WRITE(ID,99) (CEX(I),I=1,5)
C...  IF(ISW.NE.0.AND.JW.EQ.JKL) WRITE(ID,99) (CC(I),I=1,5)
C...  IF(ISW.NE.0.AND.JW.EQ.JKL) WRITE(ID,99) (AWORK(J,15),J=1,IMAXP1)
C...  STORE THE MATRIX FOR TESTING OF SOLUTIONS
C
       DO 245 I=1,IMAXP1
       DO 245 J=1,IMAXP1
 245   ASTOR(I,J)=AWORK(I,J)
C
C... CALL SIMUL TO SOLVE THE MATRIX EQUATION.
      CALL SIMUL(IMAX,AWORK,BWORK,1.0E-8,0,51,ID)
C
C...  IF(ISW.NE.0) WRITE(ID,98)
C98    FORMAT(/2X,'SOLUTIONS'/)
C...       IF(ISW.EQ.1) WRITE(ID,299) INT(ZR),(BWORK(J),J=1,IMAX)
C
C........ calculate the vibrational temperature
      CALL TVIBP(J,1.0D3,TVBN2P,BWORK)
      ITV=ITV+1
C....      IF(ITV.LT.2.AND.ABS((TVBN2P-TVSAV)/TVSAV).GT.0.1.AND.
C....     > TVSAV.GT.0.0) GO TO 23
C....      WRITE(6,*) ITV,INT(ZR)
C
      N2TOT=0.0D0
      DO 251 I=1,IMAX
 251  N2TOT=N2TOT+BWORK(I)
C........  EFFECTIVE N2+ + O->N2 + O+ RATE. PRN2P= N2* PROD BY N2+(V) +N2
      R42=RTS(42)*RTS(10)
      RAV=0.0D0
      DO 252 IR=1,5
C     PRN2P(IR,JP)=BWORK(IR)*N2*KQ
      RTS(90+IR)=BWORK(IR)*N2*KQ
 252  RAV=RAV+BWORK(IR)*RL(IR)
      RTS(99)=RAV/N2TOT
       IF(IABS(ISW).EQ.1) WRITE(ID,299) NINT(ZR),(BWORK(J),J=1,IMAX)
     > ,TVBN2P
C
C....... ALTERNATIVE CALCULATION OF N2+ DENSITY
      PROD=QX2S+QA2P+QB2S+PNOMET
      ALOSS=O*(RTS(10)+RTS(99))+RTS(11)*NE*(TVBN2P/300.)**0.39
     >  +RTS(17)*O2
      CHN2P=PROD/ALOSS
      PR(1)=QX2S+QA2P+QB2S
      PR(2)=PNOMET
      ALR(1)=RTS(10)*CHN2P*O
      ALR(2)=RTS(42)*RTS(10)*CHN2P*O
      ALR(3)=RTS(11)*CHN2P*NE*(TVBN2P/300.)**0.39
      ALR(4)=RTS(17)*CHN2P*O2
C.....      IF(IABS(ISW).EQ.1.AND.ID.EQ.3) WRITE(14,96)
C 96  FORMAT(/2X,'ALT    QX2S',4X,'QA2P',4X,'QB2S',3X,'O+(2D)',3X,'N2*+O
C    >  N2*+O  N2*+e  N2*+O2'/)
C....      IF(IABS(ISW).EQ.1.AND.ID.EQ.3) WRITE(14,299)
C....     >  INT(ZR),QX2S,QA2P,QB2S,PNOMET
C....     > ,(ALR(I),I=1,4),R42,RTS(99)
C
C... TESTING RESULTS TO SEE IF THEY ARE SOLUTIONS
C
      NOSOL=0.0
      BWORK(IMAXP1)=-1.0
      DO 255 I=1,IMAX
      DO 256 J=1,IMAXP1
 256  ASOL(I)=ASOL(I)+ASTOR(I,J)*BWORK(J)
      IF(ABS(ASOL(I)).GT.1.0E-10) NOSOL=1
 255  CONTINUE
C....      IF(NOSOL.EQ.1) WRITE(ID,97) JP,(ASOL(I),I=1,IMAXP1)
C97   FORMAT(2X,'WARNING!!   INVALID SOLUTION!'/,I6,1P,20E8.1)
C
       DO 260 I=1,4
       EMRAT(I)=A(I,1,3)*BWORK(I+5)
       IF(I.GT.2)GO TO 260
       EMRAT(I+4)=A(1,I,4)*BWORK(11)
       EMRAT(I+6)=A(2,I+1,4)*BWORK(12)
  260  CONTINUE
C 95  FORMAT(/2X,'VOLUME EMISSION RATES'/,2X,'11088',5X,'9183',5X,'NM20',
C    >3X,'NM31',3X,'3914',4X,'4278',4X,'3882',4X,'4234'/)
C...  IF(ISW.NE.0.AND.JW.EQ.JKL) WRITE(ID,95)
C...  IF(ISW.NE.0.AND.JW.EQ.JKL) WRITE(ID,99) (EMRAT(I),I=1,8)
C
 90    FORMAT(2X,1P,20E9.2)
C99    FORMAT(1X,1P,22E8.1)
 299   FORMAT(I5,1P,22E8.1)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SIMUL(N,A,X,EPS,INDIC,NRC,ID)
C........ THE ORIGINAL OF THIS PROGRAM IS IN A BOOK ENTITLED "APPLIED
C...... NUMERICAL METHODS" BY BRICE CARNAHAN, H.A. LUTHER AND JAMES
C...... O. WILKES PAGE 290. IT SOLVES A SET OF LINEAR EQUATIONS BY THE
C...... GAUSS JORDAN ELIMINATION METHOD WITH THE MAXIMUM PIVOT STRATEGY
      IMPLICIT REAL (A-H,O-Z)
      DOUBLE PRECISION A,PIVOT,AIJCK,X
      DIMENSION IROW(50),JCOL(50),A(51,51),X(51)
      DATA IROW/50*0/,JCOL/50*0/
C
      MAX=N
      IF(INDIC.GE.0) MAX=N+1
C
      IF(N.LE.50) GO TO 5
      WRITE(ID,200)
      RETURN
C
 5    CONTINUE
C5    DETER=1.
      DO 18 K=1,N
      KM1=K-1
C
      PIVOT=0.0D0
      DO 11 I=1,N
      DO 11 J=1,N
C
      IF(K.EQ.1) GO TO 9
      DO 8 ISCAN=1,KM1
      DO 8 JSCAN=1,KM1
      IF(I.EQ.IROW(ISCAN) .OR. J.EQ.JCOL(JSCAN)) GO TO 11
 8    CONTINUE
 9    IF(ABS(A(I,J)).LE.ABS(PIVOT)) GO TO 11
      PIVOT=A(I,J)
      IROW(K)=I
      JCOL(K)=J
 11   CONTINUE
C
      IF(ABS(PIVOT).GT.EPS) GO TO 13
      RETURN
C
 13   IROWK=IROW(K)
      JCOLK=JCOL(K)
C
      DO 14 J=1,MAX
 14   A(IROWK,J)=A(IROWK,J)/PIVOT
C
      A(IROWK,JCOLK)=1.0/PIVOT
      DO 18 I=1,N
      AIJCK=A(I,JCOLK)
      IF(I.EQ.IROWK) GO TO 18
      A(I,JCOLK)=-AIJCK/PIVOT
      DO 17 J=1,MAX
  17  IF(J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J)
 18   CONTINUE
C
C
      DO 20 I=1,N
      IROWI=IROW(I)
      JCOLI=JCOL(I)
C     JORD(IROWI)=JCOLI
 20   IF(INDIC.GE.0) X(JCOLI)=A(IROWI,MAX)
      RETURN
 200  FORMAT(2X,'ON TOO BIG')
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE TVIBP(J,TGUESS,TV,VPOP)
C....... to find the vibrational temperature of the equivalent Boltzmann
C....... N2* densities. Add up the energy in the calculated vibrational
C....... distribution then calculate the vibrational temperature of the nal
C....... Boltzmann distribution that has the same energy.               rgy.
C....... A Newton procedure is employed to make E(Bolt)=E(cal). TGUESS= neutral
C....... temperature for first guess, TV=calculated vibrational temperature,
C....... VPOP(51) contains the calculated N2+(V=0->5). Written by P. Richards
      DOUBLE PRECISION VPOP,ECALC,FUN,FEX,DEX,TV,H,TGUESS
      DIMENSION VPOP(51)
C........ total energy in calculated distribution = D1 E1 + D2 E2 +......
C........ E2=2*E1, E3=3*E1 etc. E1 also appears in Bolt. and cancels out
      ECALC=VPOP(2)+2*VPOP(3)+3*VPOP(4)+4*VPOP(5)+5*VPOP(6)
C........ first guess=TGUESS, H= change in TV used to calculate deerivative
      TV=TGUESS
      JITER=0
      H=0.0001*TV
C........ find function value for TV
 20   CALL VFUN(TV,VPOP(1),ECALC,FUN)
      FEX=FUN
      TV=TV+H
C........ find function value for TV=TV+H for derivative DEX
      CALL VFUN(TV,VPOP(1),ECALC,FUN)
      DEX=(FUN-FEX)/H
      IF(DEX.EQ.0) GO TO 30
C......... increment TV
      TV=TV-H-FEX/DEX
      JITER=JITER+1
      IF(JITER.GT.22) GO TO 30
C......... if DEX small enough, stop.
      IF(DABS(FEX/(DEX*TV)).GT.1.0E-3) GO TO 20
      RETURN
C......... something has gone wrong !!
 30   CONTINUE
      TV=0.0D0
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE VFUN(TV,D1,ECALC,FUN)
C......... Calculate the energy in the Boltzmann distribution temp=TV.
C......... ECALC=energy in calculated distribution, D1=density of v=0
C......... FUN= E(Bolt)-E(Cal)
      DOUBLE PRECISION ECALC,D1,EBOLT,FUN,X,X2,X3,X4,X5,TV
      X=DEXP(-2900./TV)
      X2=X*X
      X3=X2*X
      X4=X3*X
      X5=X4*X
      EBOLT=D1*(X+2*X2+3*X3+4*X4+5*X5)/(1.0-X)
      FUN=EBOLT-ECALC
      RETURN
      END
